Topological phases in ultracold polar-molecule quantum magnets 
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We show how to use polar molecules in an optical lattice to engineer quantum spin models with 
arbitrary spin S > 1/2 and with interactions featuring a direction-dependent spin anisotropy. This 
is achieved by encoding the effective spin degrees of freedom in microwave-dressed rotational states 
of the molecules and by coupling the spins through dipolar interactions. We demonstrate how one 
of the experimentally most accessible anisotropics stabilizes symmetry protected topological phases 
in spin ladders. Using the numerically exact density matrix renormalization group method, we find 
that these interacting phases - previously studied only in the nearest-neighbor case - survive in the 
presence of long-range dipolar interactions. We also show how to use our approach to realize the 
bilinear-biquadratic spin-1 and the Kitaev honeycomb models. Experimental detection schemes and 
imperfections are discussed. 
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Recent advances in ultracold polar molecules [1-3] , Ry- 
dberg atoms [4, 5], magnetic atoms [6, 7], and magnetic 
defects in solids [8-10] have spurred tremendous inter- 
est in exotic strongly-correlated many-body phenomena 
arising from anisotropic, long-ranged dipole-dipole inter- 
actions [11-43]. The types of anisotropies realizable with 
these interactions are typically limited to simple changes 
of the interaction sign and magnitude according to the 
spherical harmonic l2,o oc 1 — 3 cos^ 9, where {6, 4>) are 
the spherical coordinates of the vector connecting the two 
interacting dipoles [11, 13, 32, 33]. 

In this Letter we show, in the context of polar 
molecules, that microwave dressing provides a tremen- 
dous degree of simultaneous control over five independent 
dipole-dipole interaction terms whose angular depen- 
dences are given by the rank-2 spherical harmonics. This 
opens the door to simulating well-known models includ- 
ing the spin-1/2 XXZ model with a direction-dependent 
spin anisotropy, the spin-1 bilinear-biquadratic model 
[44], and the Kitaev honeycomb model [45]. Thanks to 
the use of direct dipole-dipole coupling, the resulting in- 
teractions are stronger and hence easier to observe exper- 
imentally than other - potentially direction-dependent - 
spin-spin interactions such as superexchange in ultracold 
atoms [46] or pcrturbative dipole-dipole-mediated cou- 
plings between polar molecules [16, 17]. 

As a specific example demonstrating the reach of our 
method, we show how to design a spin-1/2 XXZ model 
with direction-dependent spin anisotropy using a min- 
imal and experimentally reasonable microwave configu- 
ration. In a two-legged ladder geometry with nearest- 
neighbor interactions, this model has been shown to ex- 
hibit symmetry protected topological (SPT) phases [47]. 
These phases are exotic gapped states of matter distinct 
from trivial gapped phases when specific symmetries are 
present. They have recently attracted extensive interest 




FIG. 1. (color online). A lattice of polar molecules in the XY 
plane is subjected to a DC electric field along z. We define the 
xyz coordinate system as the rotation of the XY Z coordinate 
system around Z by $o and then around y by Oq. A vector R 
with polar coordinates [R, $) in the XY plane has spherical 
coordinates (PL, 6, <j)) in the xyz coordinate system. 



[48-56] because they do not fit within the framework of 
Landau symmetry breaking and possess exotic properties 
such as topologically protected edge states [57] , nonlocal 
order parameters [58, 59], and unique entanglement prop- 
erties [54, 60]. Using the density matrix renormalization 
group method (DMRG) [61], we compute the phase dia- 
gram of the two-legged-ladder model obtained in our po- 
lar molecule implementation and provide evidence that - 
at least in this one-dimensional model - SPT phases also 
exist in the presence of long-range dipolar interactions. 

In a major advance over Refs. [32, 33, 62], our proposal 
realizes an interacting topological phase. Furthermore, 
relying on homogeneous microwave - not optical - dress- 
ing, our proposal is much easier to realize experimentally 
than the one on topological flatbands [62]. Finally, since 
the relevant motional energy scale in our setup is the 
lattice bandgap, our proposal can be realized at much 
higher motional temperature than the t- J-type model of 
Refs. [32, 33], which relies on tunneling. 

Setup. — We consider an array of polar molecules con- 
fined to the XY plane and pinned in a deep optical 
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lattice with one molecule per site [see Fig. 1(a)]. Each 
molecule is treated as a rigid rotor with dipole moment 
operator d, angular momentum operator N, and a rota- 
tional constant B and is described by the Hamiltonian 
Hq — BN^ — Ed^ in the presence of a DC electric field E 
along z. As one turns on E, the simultaneous eigenstates 
of and with eigenvalues N{N -\- 1) and M adia- 
batically connect to eigenstates of Hq^ which we denote 
\N,M). The dipole-dipole interaction between molecules 
i and j separated by R [see Fig. 1(a)] is [63] 

= E {~lYCl^{0,<l^)Tl{d,,A,), (1) 

q=-2 

where Cl{6,(j)) = ^47r/5 Y2,q{e,(j}) and the many-body 
Hamiltonian is iJ = (1/2) E.^^ Here T'^{di,dj) 

is given by = dfdf, r|i = {d'^df + dfd"^) /V2, and 
T2 = (drd+ + 2d'ld° + d+d-j) /^/6, where rf° = d' and 

d± = =F(d'' ± idy)l^/2. Thus changes the total M of 
the two molecules by q. For R ^ 0.4 fim, the interaction 
energy scale is d'^/R^ - 1 (100) kHz in KRb (LiCs). 

To obtain a spin-S Hamiltonian, we select in each 
molecule 2S'-|- 1 disjoint sets of jiV, M) states and couple 
the states within each set to form dressed states (^ n 
microwave fields are needed to couple n states) . We then 
choose one [64] dressed state from each set to create the 
spin-5 configuration. Projecting Eq. (1) onto the cho- 
sen spin-S* basis, the resulting spin-spin interactions con- 
sist of five potentially independently controllable terms 
with angular dependences Cq, Re[Ci], Im[Cf], Re[C|], 
and Im[C|]. Refs. [32, 33] considered the special case 
oi S = 1/2 in the presence of the Cq term alone and 
limited the discussion to situations where total was 
conserved. In this Letter, we evince the power of the 
approach beyond this special case. 

Interactions featuring a direction- dependent spin 
anisotropy. — Our first demonstration of novel direction- 
dependent interactions focuses on 5 = 1/2 and as- 
sumes that Hij connects a pair of molecules in the state 
|mi)|m2) (jn^i) and |m2) are dressed states) only to it- 
self and to |m2)|mi), while all the other processes are 
off-resonant and thus negligible. Although one may be 
able to independently control each of the five terms, 
here we will focus on the Cq_|_2 terms since C^.^ terms 
are resonant only at specific values of E [65] . 

Consider the level configuration in Fig. 2(a). We as- 
sume ri± > and |ri±| > Hij and take jf) = |0, 0) and 
11) = ajl, —1) — /3|1, 1) as our dressed spin states, where 

{a,/3} = {n_,n+} /^nl, +nl. Notice that \i) is a 

single-molecule eigenstate in the presence of fl±. The 
spin model is then derived by projecting Hij onto states 
It) and \l) via the same steps as in Ref. [33] with one 
major difference: dfdj', featuring a C^2 angular depen- 
dence, resonantly couples |1,-1)|0,0) |0,0)|1,1). The 
result is [66] 

R'H,, = M'^)S^S] + J,yi-f){S:S^ + SfS^), (2) 




FIG. 2. (color online), (a) The level scheme and reso- 
nant microwave coupling used to realize the Hamiltonian, 
Eq. (2). The dressed states we choose are {|t) , |i)} = 

{|0,0) , (n_ |1,-1) - fl+ |1, 1))/^QPl +02.}. (b) Microwave- 
dressed rotational levels for the SU(2)-symmetric spin-1 
model. The dressed states are 1 1) (linear combination of states 
indicated by triangles), |0) (ovals), and | — 1) (the rest). The 
diagram is schematic: the real system is anharmonic and lev- 
els \N,M) with the same A'' are non-degenerate (unless the 
levels have the same \M\). 



where Jz{^) = (1 - 3cos2(* - 4>o)sin^ 6o)(mo - Mi)^: 
J,y($) = -/i2^(l - 3cos2($ - $o) sin^ Oo) + &aP^ll^[l - 
cos2($ - $o)(1 + cos^eo)], Mo = (0,0|dO|0,0), ^ll = 
(l,l|d°|l,l), noi = (l,l|d+|0,0), and Sj is the spin-1/2 
operator for molecule j. The spin anisotropy Jxy/Jz of 
this XXZ model changes depending on the polar angle $ 
of the vector R connecting the two interacting molecules. 
As we discuss below, Eq. (2) allows one to study SPT 
phases in ladders. Another special case is a square-lattice 
Heisenberg model with a tunable ratio between coupling 
strengths on X and Y bonds [67] . In the nearest-neighbor 
limit, this enables one to study the change from one- 
dimensional (uncoupled) chains to a two-dimensional be- 
havior. Such models have also been used to explore 
the physics of stripes in high-temperature superconduc- 
tors [68, 69]. While we see that even the simple level 
structure of Fig. 2(a) yields a wealth of exotic physics, 
we show in the next two sections that additional features 
can be accessed with increased microwave control. 

Degenerate dressed states and non-Abelian anyons. — 
To realize models such as the quantum compass model 
[70], the Kitaev honeycomb model [45], and the Yao- 
Kivelson model [71], we need to go beyond Eq. (2) and 
realize terms, such as SfSj, that do not conserve the 
total S^. To do this, we simply tune |t) and jt) to be 
degenerate. 

As an example, consider the Kitaev honeycomb 
model, where interactions along $ = 7r/6, 7r/2, and 
57r/6 are of the form SfSJ, Sf , and SfS^, re- 
spectively [45]. At (Oc'l'o) = (0,0), the interaction 
between two molecules i and j is R'^Hij{^) ~ v(<I>) • M, 
where v($) = {1, -3cos(2$)/2, 3sin(2$)/2} and 
M = {^Tl Ti + Tl^,zTi - zT^J. Since v(7r/6), 
v(7r/2), and v(57r/6) are linearly independent, it is, 
in principle, possible to choose the degenerate dressed 
states It) and |t) to ensure that Hij{TT/6) oc SfSj, 
H,j{TT/2) cx SfS], and H,j{5n/6) cx S^S^. In Ref. [72], 
we show that - with sufficient microwave control - such 
a choice of dressed states is indeed possible, allowing 
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one to realize the Kitaev B phase in the presence of a 
magnetic field. This gapped phase supports non-Abelian 
anyonic excitations, which can be used, for example, for 
topologically protected quantum state transfer [73] and 
quantum computing [45]. 

> 1/2 and the bilinear-biquadratic model. — We now 
show that one can extend this tremendous control over 
spin-spin interactions to S" > 1/2. In particular, we 
show how to obtain the general SU(2)-symmetric spin- 
1 Hamiltonian, i.e. the bilinear-biquadratic Hamiltonian 
cos(7)Si • Sj + sin(7)(Si • Sj)^, which has a rich phase 
diagram even in one dimension [17, 44, 74]. In particu- 
lar, 7 = 7r/4 and arctan(l/3) give the SU(3)-symmetric 
and the AKLT [75] Hamiltonians, respectively. One can 
also consider generalizations to SU(iV) with arbitrary N 
as well as away from SU(2) [76]. 

We build our spin-1 dressed-state basis {|1) , |0) , | — 1)} 
from the 15 bare levels shown in Fig. 2(b). For simplicity, 
the model presented here will have only the Cg term. We 
work at E = 3.2UB/d [= 13 (7) kV/cm in KRb (LiCs)], 
for which the bare-states process |1,0) |1, 0)— >-|2, 0) |0, 0) 
is resonant. Furthermore, we choose the energies of the 
dressed states to make the process |0) |0) | — 1) |1) reso- 
nant [77]. The latter is needed to engineer the Sj' 
term present in the desired Hamiltonian. Aside from 
this exception, we again assume that a pair of molecules 
in dressed states \mi)\m2) is connected via Hij only to 
itself and to \m2)\mi). Using 12 microwaves to create 
the dressed states in Fig. 2(b), we find [66] that we can 
achieve any 7 and, thus, any bilinear-biquadratic Hamil- 
tonian. Removing all five = 4 states, we are left with 
just 8 microwaves, which simplifies the experimental im- 
plementation but at the cost of only accessing the 7 = 1.1 
point. The typical strength of interactions achieved is 
R^Hij ^ O.Old^. Stronger interactions and a reduced 
number of microwaves might be achievable by further 
optimizing the choice of levels and microwaves. 

SPT phases in spin ladders. — We now turn to the fo- 
cus of our work: we use Eq. (2) to implement a spe- 
cific ladder model [see Fig. 3(a)] introduced in Ref. [47] 
and shown to support nontrivial SPT phases for nearest- 
neighbor interactions. The symmetries protecting these 
interacting topological phases are the exchange a of the 
two legs and D2 = {E, Rx,Ry,Rz}, where E is the iden- 
tity and Ra is a 7r-pulse around the axis a on all spins. 
Note that although we focus on a ladder system since it 
is amenable to a numerically exact treatment, we expect 
an even richer phase diagram in dimensions D > 1. 

While for our choice of levels b = — /^i)^ sat- 

isfies b e [2.6, cx)), any 6 > can be accessed by us- 
ing reduced nuclear spin overlaps between [f) and |^) 
(see "Experimental considerations" below) or by choos- 
ing \N,M) with different A^. To ensure the cr symme- 
try, we take — tt/2. To ensure a Heisenberg Hamil- 
tonian along $ = 0, we choose a/3 = {b + l)/(65). 
Since a/3 e [0,1/2] and b > 0, b can go from 1/2 
to 00. Rescaling the interaction by (/ig — Mi)^ [which 
goes up to w O.ld^ for our choice of levels], defining 



A2 = 1 — 3cos^ Qq (tunable via Gq between —2 and 1), 
and \xy = _ ^(''+i) „ 4b+i ^^ [see shaded area in Fig. 3(b)], 
we obtain Eq. (2) with Jz($) = 1 - (1 - A^) sin^ $ and 
Jxyi^) ~ 1 — (1 — \xy)sm^^. As desired, at nearest- 
neighbor level, these expressions for Jz{^) and Jxy{^) 
reproduce Fig. 3(a). 

In our implementation using polar molecules, the 
nearest-neighbor interactions are replaced by dipolar in- 
teractions, which give rise to nontrivial longer range cor- 
rections. In order to investigate the role of these cor- 
rections, we numerically calculate the phase diagram of 
the spin ladder with long-range interactions, as shown in 
Fig. 3(b), using DMRG on 200 rungs with smooth bound- 
ary conditions [78] . By performing a finite-size scaling us- 
ing systems with up to 400 rungs, we estimate the finite- 
size effects to be comparable to the size of the symbols in 
Fig. 3(b). The phase diagram is qualitatively similar to 
the nearest-neighbor case [47] and exhibits the same six 
phases, including the two SPT phases. In Ref. [47] 's Ian- 
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FIG. 3. (color online). Spin ladder of Ref. [47] and its phase 
diagram showing SPT phases in the presence of long-range 
interactions, (a) The nearest-neighbor model, (b) Phase dia- 
gram in the presence of long-range interactions. The shaded 
area indicates points achievable with the simple configuration 
of Fig. 2(a). The shaded area does not extend past the lim- 
its of the vertical axis, while it extends infinitely far along 
the horizontal axis, (c) Entanglement splitting (open boxes, 
left axis) and energy gap (solid boxes, right axis) of the SPT 
phases along = —0.5 cuts (bold, red lines) shown in (b). 
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guage, the four nontopological phases are the Ising Ncel, 
the Ising stripe Neel, and two product phases of rung 
singlets and 5^ = rung triplets. The remaining two 
phases, to and t^, are two out of seven nontrivial SPT 
phases protected hy D2 x a [47]. The to phase can be 
connected to the Haldane phase [79] and to the AKLT 
state [75] by treating the triplet states on each rung as 
a spin-1 particle. Meanwhile, in the nearest-neighbor 
case, is obtained from to by taking )t) — |t) on 
one of the legs, which is the Xxy — >■ —Xxy symmetry of 
the nearest-neighbor phase diagram [47]. Long-range in- 
teractions break this symmetry and, in particular, reduce 
the size of the phase relative to the to phase as a re- 
sult of substantial ncxt-nearest-neighbor SfSJ + SfSj 
interactions for X.j.y > 0. 

We observe that the t^ phase is sensitive to artificial 
cutoffs in the interaction range, so we use matrix product 
operators within DMRG to provide an efficient descrip- 
tion of interactions without a cutoff: instead, we fit the 
long-range interactions to a sum of exponentials [66, BO- 
SS]. The boundaries of the Ising Neel and Ising stripe 
Neel phases were obtained by calculating the correspond- 
ing order parameters. We identified the two rung phases 
by calculating {SfSJ + SfS^^) on the rungs and by ver- 
ifying that the gap doesn't close as \\xy\ increases to 
large values where the system is ultimately exactly solv- 
able. The boundaries of the SPT phases were obtained by 
computing the entanglement splitting, which we define as 
ES = X]j=odd (^i ^ where Wj are the eigenvalues 

of the reduced density matrix for a bipartition at the cen- 
ter of the system, sorted from largest to smallest. Due 
to their two- fold degenerate entanglement spectrum [60] , 
SPT phases have ES — 0, as shown in Fig. 3(c). We have 
also verified that all phases are gapped. Interestingly, 
the energy gap in the SPT phases, shown in Fig. 3(c), 
exhibits a cusp indicative of a level crossing, which de- 
serves further investigation. Finally, using again systems 
with 200 rungs, we added to the Hamiltonian a small 
term oc Sf ± 5J, where i and j are sites on the same edge 
rung. Referring to operators that split the edge degen- 
eracy as active operators, we verified the prediction [47] 
that Sf + Sj is an active operator for the to phase while 
Sf — S'j is not, and vice versa for the tz phase. 

Experimental considerations. — As suggested in 
Ref. [47], an SPT phase can be classified by finding 
its active operators, i.e. those operators that split the 
edge degeneracy. We propose to diagnose this splitting 
by measuring an active operator in linear response to 
the application of that same operator at frequency w 
and looking for the zero-bias (w = 0) peak. Repeating 
the same procedure for inactive operators will yield no 
zero-bias peak. In our implementation, a z magnetic 
field proportional to /Zg — /xf naturally arises at the edges 
from dipole-dipole interactions [66]. Since such a field 
constitutes an active operator of the to and tz phases, it 
is natural to probe the response of the system by tuning 
/ig — with the DC electric field. In combination 
with a spectroscopic verification of the bulk gap, the 



response to active operators allows one to detect and 
classify SPT phases. A more modest first experimental 
step could be to use a Ramsey-type experiment [40] to 
benchmark how accurately the molecules emulate the 
desired Hamiltonians. 

Polar alkali dimers have hyperfine structure H^f [33, 
84], which we have ignored so far. We will illustrate 
how to deal with H^f for the specific case of /S = 1/2. 
Assuming that microwave Rabi frequencies fij are much 
larger than H^f, we can project the hyperfine structure 
on dressed states Jf) and \\.). The necessary conditions 
i?hf <C ^ i? are easy to satisfy: for example, in 
"OK^^Rb, Hhf ~ (27r)l MHz and the rotational constant 
is B ~ (27r)l GHz. The simplest situation arises when 
an applied magnetic field - of a strength already ex- 
perimentally used [1] makes (t| i?hf |t) and (|) Hm ]|) 
diagonal in the same basis of decoupled nuclear spins. 
The nuclear-spin degree of freedom can then be elimi- 
nated by working with a single state from this basis. For 
smaller magnetic fields, one could prepare the system in 
any pair of non-orthogonal eigenstates of (t| i^hf It) and 
{l\ Hhf II). An imperfect overlap of these two states will 
effectively reduce the transition dipole moment between 
It) and I4-), resulting in an additional control knob of the 
interactions. 

Controlling tens of independent microwave frequen- 
cies in the frequency range required by our proposal is 
straightforward [85]. The two uncertainties involved are 
in the generation of the microwaves and in the coupling 
to molecules. The latter is dominant: current ultracold 
molecule experiments observe only 0.1% deviations in 
their ~ 1ms microwave pulses without any particular op- 
timization [86] , and this is expected to be independent of 
the number of microwaves applied. Polarization control 
is more challenging. However, this should be attainable, 
for example, by simply interfering the outputs of two in- 
dependently controlled microwave horns. 

Outlook. — While dipolar interactions did not destroy 
the SPT phases in our example, quantum magnets with 
long-range interactions have recently been shown to har- 
bor novel, unusual, and often dimension-specific physics 
[87-93]. The polar- molecule experiment we propose 
could therefore help guide the theoretical understanding 
of these effects in 2D and 3D systems - including SPT 
phases - where efficient numerical methods are not avail- 
able. In fact, the classification of SPT phases is yet to be 
extended to models with long-range interactions. Finally, 
we expect our methods to be immediately extendable to 
other dipole-dipole interacting systems such as Rydberg 
atoms [4, 5], magnetic atoms [6, 7], and magnetic defects 
in sofids [8-10]. 
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I. DETAILS ON THE MOLECULAR PHYSICS 

A. Derivation of Eq. (2) in the main text 

Defining {|0), |1), |2)} = {|0, 0), |1, -1), |1, 1)} and 
nm = the dipole-dipolc interaction between 

molecules i and j that are R = [R, 9, <f)) apart is 



R^Hij = (1 - Scos^e) 



(SI) 



{^lono + Mini + Min-2)(Mo«o + Mi^-i + ^1^2) 



-i/i2i(|01)(10| + |02>(20|+h.c.) 

+Mgi|sin2 [e^2'^(|01)(20| + |10)(02|) + h.c] . 

Projecting on |t) = |0) and ||) = a | — 1) — keeping 
only the terms that conserve the total , and dropping 
a constant, we obtain Eq. (2) in the main text with an 
additional term (l-3cos^($-$o) sin^ Qo){^J'o - IJ-Di^t + 
Sj)/3 on the right-hand side. For a lattice with a single- 
site or two-site (as in the ladder) unit cell, neglecting 
edge effects, the additional term is a uniform magnetic 
field, which is irrelevant since J2i'^i conserved [SI]. 
On the other hand, if one is interested in studying edges 
in the presence of nonzero Jz, one can choose a level 
configuration where /io = —Mi 7^ 0, so that the addi- 
tional term vanishes. This is achieved, for example, for 
{|t),|i)} = {|l,0),a|2,~l)-/3|2,l)} at dE/B = 5.072 
or for {It), 14)} = {|2,0),a|2,-l) -/3|2,1)} at dE/B = 
10.535. Since a z magnetic field is an active operator for 
both the to and the phases [S2], tuning E away from 
the value where Mo = —Mi allows one to controllably 
turn on active operators and, thus, probe the nature of 
the edge states, as discussed in the main text. 



B. Details behind Fig. 2(b) in the main text 

We number the states in Fig. 2(b) in the main 
text in the left-to-right and bottom-to-top order as 

{|1),|2),|3),...,|15)} = {|0,0),|1,0),|1,1),...,|4,4)}. 
We then write the three dressed states as \p) = 
T.j(p) y/Xj\j), where Xj > 0, T,j(p)Xj = I, p = 0,±1, 
and j{p) means that j runs only over the 5 states be- 
longing to \p) in Fig. 2(b). Whether |1) refers to the 



bare state or to the dressed state will be clear from the 

context. Four microwave fields coupling five bare states 
that make up each dressed state allow one to arbitrarily 
tune the composition Xj and energy of the dressed state. 
Thus 12 microwaves are needed to fully control all If) Xj. 
Keeping only resonant terms [as in Eq. (SI)] and pro- 
jecting onto dressed states, the dipole-dipole interaction 
between two molecules that are R = {R, 9, (j)) apart is 



1-3 cos^ ( 



i?3 



^ ApA^ \pq) {pq\+J2 Bp\PP) (PP\ 



L p,g 



E +11.C.) + J+(|00)(-11| +h.c.) 



p<q 



.(S2) 



where Ap — J2j(p) -^iMji -^p — J2i(p)<j{p) XiXjdij, Jp^q 



dij 



2m| ifMii)=M{j) 
-4 if\M{i)-M{j)\ = l , (S3) 
otherwise 



= and M^ = {i\d°\i). To allow for J+ 

to be negative, we can simply change the sign of ^/x4 in 
the expansion of dressed state |1). 
To achieve 

R^Hij = (1 - 3cos^ 9) [ai + aaS, • + 03(8, • Sj)^] (S4) 

we find the constraints on Ap, Bp, Jp^q, and J+ in terms 
of Qi by matching the matrix elements in Eqs. (S2) 
and (S4). It is straightforward to numerically check 
that, by tuning xj, these constraints can be satisfied 
for arbitrary 7 = arg(a2 + ia^). As an example, setting 
Xj = for J = 11,..., 15, we are left with just 8 
microwaves (two microwaves are needed to couple |2, 1) 
to |3, 3)), and they give us a solution with {02,03} = 
{0.0027,0.0055} (i.e. 7 = 1.1) for {xi, ^12, . . . , Xio} = 
{0.0189, 0.1575, 0.2342, 0.5477, 0.9654, 0.7132, 0.0209, 
0.1293,0.1972,0.0157}. 

It is worth pointing out that, in many geometries, 
5'f + and (5'f )2 + (S'|)2 in Eq. (S4) just lead to uni- 
form single-particle energy shifts in the bulk. Such shifts 
can be easily compensated by tuning the energies of the 
dressed states. Thus, in cases where one is not worried 
about introducing 5? and {S^)'^ terms near the boundary. 
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the conditions used above can be relaxed even further. 
For example, we can then set Xj = for j — 9, . . . , 15, 
and we find a solution requiring only 5 microwaves - 
with {02,03} = {0.0344,0.0041} (i.e. 7 = 0.12). 



II. DETAILS ON THE NUMERICS 

All density matrix renormalization group (DMRG) cal- 
culations were performed on finite ladders with open 
boundary conditions in both the X and Y directions. 
Our calculations conserved total S^, and ground states 
were found by targeting the S"^ = sector. In order to 
attain discarded weights less than 10~^, we typically kept 
between 100 and 1000 density-matrix eigenstates. 

To work with long-range interactions in DMRG, we 
used the ITensor library [S3] and expressed our Hamilto- 
nians as matrix-product operators, which can exactly en- 
code long-range exponentially decaying interactions. We 
then approximated our dipolar interactions, which decay 
as as a sum of exponentials. (For technical details 

of this procedure, see below.) Typically, using just 5 ex- 
ponentials on systems of up to 400 nmgs was sufficient 
to obtain fits deviating from the exact interactions by 
less than 10~^ at any fixed range. We also checked con- 
vergence of observables: for example, the entanglement 
splitting of 400-rung systcmis fit with 5 exponentials dif- 
fered from results with 12 exponentials by less than 10~^. 

The data in Fig. 3(c) of the main text was obtained 
via finite-system calculations on system sizes typically up 
to 400 rungs, extrapolated to the thermodynamic limit 
using quadratic fits. All extrapolations were of good 
quality (errors order of symbol sizes) except for the en- 
tanglement splitting at X^y = 1.04, which extrapolated 
to a negative value. The reported value is half of the 
value of the largest system. To calculate energy gaps 
for Fig. 3(c), we computed the lowest-lying = and 
= 1 states, sequentially orthogonalizing against lower 
lying states within DMRG [S4]. 

The data in Fig. 3(b) was primarily obtained on fixed- 
size systems of 200 rungs. However, to reduce finite- 
size effects, we employed smooth boundary conditions 
[S5], smoothly reducing all spin-spin interactions from 
full strength to zero over a region of ss 20 rungs at each 
edge. In this way we were able to obtain results similar 
to open-boundary systems of roughly twice the size. 

To verify that we correctly identified the rung phases, 
we calculated the correlation function Cxy = {SfSj + 
SfSj) on the rungs. As expected, we found that Cxy < 
{Cxy > 0) in the rung singlet (triplet) phase, and \ Cxy \ — >■ 
1/2 in both phases as \Xxy \ —t- 00 at a fixed A^. 

The Ising Neel and Ising stripe Neel phases with well- 
defined order parameters were obtained without intro- 
ducing a symmetry-breaking field. Instead, the choice of 
the initial state for the DMRG algorithm picked out one 
of the two ground states, while the other one could be 
obtained by flipping all the spins in the initial state. 



A. Exponentially Decaying Long-Range 
Interactions with Matrix Product Operators 

In conventional implementations of the DMRG algo- 
rithm, the Hamiltonian is treated as a sum of individual 
terms with each term projected separately into the local 
basis used for each DMRG step [S6]. If the Hamilto- 
nian contains long-range interactions, this approach re- 
quires including each pairwise interaction term such that 
DMRG no longer scales linearly in system size even in 
one dimensional gapped phases. 

DMRG has been understood to be a method for opti- 
mizing variational wavefunctions known as matrix prod- 
uct states (MPS) [S7, S8]. The MPS form of the wave- 
function suggests that operators may be written in a sim- 
ilar form, known as a matrix product operator (MPO): 



W = 



E 

{s,s'},{a} 



(S5) 



Besides being more convenient to use in MPS-based al- 
gorithms, rewriting the Hamiltonian as an MPO offers 
key technical advantages. Most remarkably, a finite- 
bond-dimension MPO can exactly represent a Hamilto- 
nian with exponentially decaying long-range interactions 
[S9, SIO, Sll]. 

As a concrete example, the following MPO consisting 
of only 3x3 matrices 



SI 







Jx Sj tj S j I j 



(S6) 



encodes the exponentially long-range interacting Ising 
model: 



H 



N 
i<j 



(S7) 



(Setting A = restores the conventional nearest-neighbor 

model.) Here we notate MPO tensors as matrices 
of operators. For example, Eq. (S6) indicates that 

(M^i i)^^4- = {IjY^^'i and (1^2,1)"'"^ = iS^Y''''- We also 
assume open boundary conditions, taking only the last 
row of Eq. (S6) on site 1 and the first column of Eq. (S6) 
on site N. Because the MPO has a finite, system-size- 
independent bond dimension, it can be used to study 
exponentially decaying long-range interactions within 
DMRG while retaining linear scaling in system size (in 
ID gapped phases). 



B. Power-Law Decaying Long-Range Interactions 
with Matrix Product Operators 

An MPO of finite bond dimension cannot exactly 
represent a Hamiltonian with power-law decaying long- 
range interactions. However, for a large enough expo- 
nent 7, power-law interactions can be well approximated 
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on a finite-size system as a sum of A'exp exponentials 
[SIO, Sll]: 



Power-Law Decaying Long-Range Interactions 
for Ladder Systems 



iVe: 



(S8) 



with A^cxp fairly small. We use the particularly elegant 
fitting procedure described in Ref. [Sll]. 

Again using the Ising model as an example, the Hamil- 
tonian 



N 

i<j 3 



can be represented within the approximation (S8) by the 
following MPO 




Xilj 
Aa/j 



-hS^ JxiSI Jx^si 



(SIO) 



More compactly, this can be written in block-matrix form 

as 







XSI {IX) I J 
-hS^ Jxsi Ii 



(Sll) 



noting the resemblance of the above to Eq. (S6). 



The standard approach for applying DMRG to ladder 
systems is to map the sites into one dimension through 
the following ordering: 




We may then study long-range interacting ladder systems 
within DMRG by applying and generalizing the MPO 
techniques discussed above. 

For example, to encode the same long-range interac- 
tions along each leg, but not between the legs, we allow 
the MPO matrices to be different on each leg: 



Ij 




{lX)Ij 
Jx5? 



(S12) 



w. 



(S13) 



Ij 0' 

(2) _ {lX)Ij 



JxS'l 

A different but related pattern gives only inter-leg inter- 
actions. We can include both intra- and inter-leg inter- 
actions by combining both patterns in a block-diagonal 
form. 

In this work, we consider Eq. (2) in the main text 
with J,X^) = 1 - (1 - A.^) sin^ $ and ./,.,,($) = 1 - (1 - 
Xxy) sin^ where $ is the polar angle of the vector con- 
necting the two interacting spins, as shown in the figure 
above. 

To fit these interactions to sums of exponentials, we 
must perform three fits - (Heisenberg) intra-leg, zz inter- 
Icg, and xy intcr-lcg interactions - to yield the three sets 
of fitting parameters: (x, A), {xz,K), and {Xxy,Xxy). 

With these definitions, the MPO matrices W^^'' on leg 
1 and on leg 2 take the form [S12]: 
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(S14) 
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X.A.^I 









{IXxy)Ij 



2Xxy^xySj 



{IXxy)I^ 



xyj-^j 
2Xxy^xySj Ij 



(S15) 



[SI] We will either be working at a constant total or will 
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